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Abstract 

A quite unusual diffuse scattering phenomenology was observed in the single-crystal X-ray 
diffraction pattern of cubic perovskite BMT (BaMg 1/ / 3 Ta2/303). The intensity of the scatter- 
ing is parametrized as a set of cube- like objects located at the centers of reciprocal space unit cells, 
resembling very broad and cubic-shaped (l/2,l/2,l/2)-satellites. BMT belongs to perovskites of 
formula AB'^gB",^ (A=Mg, B' =Ta, B" =Mg). The cubes of the intensity can be attributed 
to the partial correlations of the occupancies of the B site. The pair correlation function is the 
Fourier transform of the diffuse scattering intensity and the latter's idealized form yields the un- 
usual property of a power-law correlation decay with distance. Up to now this is observed only in 
a few exotic instances of magnetic order or nematic crystals. Therefore it cannot be classified as a 
short-range order phenomenon, as in most situations originating diffuse scattering. A Monte-Carlo 
search in configuration space yielded solutions that reproduce faithfully the observed diffuse scat- 
tering. Analysis of the results in terms of the electrostatic energy and the entropy point to this 
phase of BMT as a metastable state, kinetically locked, which could be the equilibrium state just 
below the melting point. 

PACS numbers: 77.84.-s Dielectric, piezoelectric, ferroelectric, and antiferroelectric materials; 61.05.C- X- 
ray diffraction and scattering; 61.43.Bn Structural modeling: serial-addition models, computer simulation; 
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I. INTRODUCTION. 



An experiment using X-ray single crystal diffraction at the Swiss-Norwegian Beam Lines 
(SNBL) of European Synchrotron Radiation Facility (ESRF) on a single crystal of cubic 
perovskite BMT (BaMg 1//3 Ta2/303) has shown an interesting diffuse scattering pattern, con- 
sisting of cube-shaped domains of nearly constant intensity, with the centers located on 
the (1/2,1/2,1/2) points of the reciprocal lattice and with an edge of ~ 1/3 of a reciprocal 
lattice unit (rlu). All the cubes have the same intensity level, apart from the usual slow 
modulations due to atomic and thermal factors. 

BMT belongs to perovskites of formula AB^gB^Og (A=Ba, B' =Ta, B" =Mg) where 
different species on the B-site can be ordered to various degrees. The scattering depends on 
the site occupancy and when the ordering is directly related to the pair correlation function 
of the different atomic sites pQ. In fact, BMT in its cubic phase is necessarily intrinsically 
disordered in the B-site. In this paper we shall mainly explore the chemical order at the 
B-site and show that the cube-like diffuse scattering can be fully explained in this way. 

Lattice deformations and displacements, due to different atomic radii of Mg and Ta 
(see Ref. [2]), are expected. So there is additional diffuse scattering originating from these 
correlated displacements. In fact streak-like diffuse scattering - relatively weaker - has been 
found at the foot of the cubic Bragg peaks, see Fig. [Tj We shall address also this point in a 
future paper. 

The simplest possibility for B-site ordering is the formation of subregions enriched in 
either Mg or Ta cations. This possibility, however, has been thoroughly ruled out as en- 
ergetically unfavourable in AB^B^O^ perovskites 0, E] because of the large charges in 
the subregions. There are then different homogeneous ways of ordering the B'-B" cations. 
These are: 

i) homogeneous random distribution of B-cations; 

ii) correlated distribution with occupation numbers being a well-defined function of dis- 
tance between two lattice nodes, but without periodicity; 

iii) periodic long range ordering of the Mg/Ta concentrations in different sublattices. 

Situation i) would give an almost wavevector-independent diffuse scattering, but this is not 
observed. B-site ordering of type iii) has been observed in BMT, with simultaneous lowering 
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of the cubic symmetry. In fact, BMT may exist in an ordered trigonal structure (P3ral) 
modification [21 [5HH] and it has been shown that this modification [2] mainly entails perfect 
periodic ordering of the B-site. The trigonal lattice is essentially a threefold superlattice of 
the cubic one, with extremely small geometric distortions, realized by stacking periodically 
a sequence of one Mg plane and two Ta planes orthogonally to one of the cubic (1, 1, 1) 
directions. The atomic sites are also very close to those of the underlying cubic sublattice, 
with only small atomic displacements consistent with the reduced symmetry. It is notewor- 
thy that the trigonal BMT superstructure - usually cast unto four equally populated domain 
orientations with the threefold axis along either of the cubic (111) directions, so that the 
average symmetry is conserved - can be considered as a commensurate distortion of cubic 
BMT. In this framework, the modulation contributes eight additional (±1/3, ±1/3, ±1/3) 
satellite Bragg peaks to the diffraction pattern, and these coincide with the vertices of the 
diffuse scattering cubes. Another partial B-site ordering type of class iii) has been ob- 
served [3111] in other AB' 1/3 B^ /3 3 cubic perovskites (PMT, with A=Pb, B' =Ta, B" =Mg). 
These have (1/2,1/2,1/2) satellites and correspond to a partial B-site ordering with a 2x2x2 
fee superstructure unit cell. The B-sites are of two types, one of pure B" and the other con- 
taining the rest of B" and all of B' cations, and the two B-site types arrange in a NaCl-type 
structure that consists of two inter-penetrating fee lattices. These satellites coincide now 
with the center of the diffuse scattering cubes observed for BMT. 

Diffuse scattering cubes straddle the reciprocal space location of the (1/2,1/2,1/2) satel- 
lites and the (±1/3, ±1/3, ±1/3) satellites identifying the two kinds of periodically ordered 
type iii) B-site structures. Although this diffuse scattering is not Bragg-like, in a sense that 



will be specified in Sec. |III B[ this suggests that the diffuse cubes result from B-site ordering 
as well, possibly from a non-periodic order of type ii). 

We stress here that the observed cubes of diffuse intensity are not ideal, as careful inspec- 
tion (c/. Fig. [T]) shows that "cubes" have a minimum in the center about 10% lower than at 
the periphery and the corresponding intensity distribution is not a flat-top. We see in par- 
ticular that intensity tends to increase near the vertices (the (±1/3, ±1/3, ±1/3) satellites 
of the trigonal superstructure) and decrease in the center - the (1/2,1/2,1/2) satellite of the 
hypothetical partly ordered NaCl-type superstructure, indicating an embryo of transition to 
the trigonal phase. We must therefore consider the phase with idealized cubes as a limiting 
situation that may not be perfectly realized in our sample. As we argue in the following that 
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the phase with ideal diffuse cubes is not the ground state at room temperature but only 
a metastable phase, kinetically hindered from transforming to the trigonal ground state, 
this is perfectly possible, as a germ of transition may have occurred due to finite quenching 
speed. However, in this framework, the limiting phase with ideal diffuse cubes is a key to 
understanding the B-site ordering phenomena in BMT, and therefore we will deal mainly 
with this phase. 

A. A few remarks 

As we are discussing the multiple phases of BMT, we should comment on an important 
point about phase stability. In the thermodynamic sense, a stable structure observed in 
an experiment corresponds to a minimum of the free energy at the conditions of crystal 
formation. At ambient conditions the Mg-Ta distribution is quenched and experimentally 
observed configurations may not necessarily correspond to a true ground state. A disordered 
or incompletely ordered state can be stabilized at high temperature by a strong entropy 
term and could be a metastable state at lower T. In fact, at room temperature the diffusion 
of Mg and Ta is practically nonexistent, so that the non-ergodic effects in the scattering 
experiments can be neglected. 

If the first Born approximation is valid and the scattering can be considered as elastic, 
the scattered intensity is the Fourier transform of the scattering density's pair (two-body) 
correlation function [9HTT] . The scattering density is the atomic structure, and for X-rays it 
is the sum of the single atom electron densities and the only degrees of freedom come from 
the spatial distribution (structure) of the atoms themselves. The pair correlation function 
(crystallographers name it the "Patterson function") is the information contained in the 
scattering pattern. A natural ambiguity arises in the relationship between the atomic struc- 
ture and the pair correlation function, as there may be many atomic structures that have the 
same pair correlation function [121 H3] • In some cases, this is not an accidental degeneracy. 
When the internal energy U depends only on the pairwise interactions, different structures 
with identical pair correlation functions are isoenergetic. At any finite temperature T, the 
multiplicity of configurations increases the entropy S and this in turn decreases the free 
energy, G = U — TS, to the advantage of the stability of the (statistically defined) structure. 
We intend to show that the B-site order in BMT, that gives rise to the cubes of the diffuse 
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intensity belongs to this case, based on the fact that pair (electrostatic) interactions between 
the B cations are sufficient to explain the phase diagram of the different BMT phases which 
derive from the ordering of B-sites. 

The pair (two-body) and the set of the higher order (many-body) correlation functions 
can fully describe any structure in an unambiguous way, including stochastically defined 
ones [U [13] . Considering higher order correlation functions is necessary, however, only when 
they make important contributions to the internal energy. This is sometimes the case, for 
instance, for structures that show displacive diffuse scattering [HHTU]. In any case, higher 
correlation functions are also a way of describing a structure in the crystal-chemical sense. 
This is especially convenient when the structure itself cannot be easily represented because 
it is only statistically defined. However, simpler ways of presenting the crystal-chemical 



informations are also possible, see Sec. HI G 



Our approach is as follows. The pair correlation function that is related to the diffuse 
cubes can be simply analyzed in terms of the B-site chemical order as a site concentration 
pair correlation function [TJ. This turns out to have a simple analytical form that shows 
a rather peculiar power-law decrease with distance that is observed in only a few cases for 
rather exotic systems. 

Monte- Carlo simulations were performed for a large scale in order to find the B-site config- 
urations that give the same pair correlation as the diffuse intensity. Successfully determined 
configurations were then analyzed in terms of a long-range order parameter showing that 
these partially ordered configurations that are compatible with the observed scattering phe- 
nomena are not distinguishable from completely disordered systems. Next, we consider that 
BMT is an ionic compound and that the electrostatic energy is the main term in the crystal 
Hamiltonian. The electrostatic energy of a rigid system depends only on pair interactions 
and the electrostatic energies and the entropy terms for the different BMT phases were cal- 
culated. These gave a good correspondence with the observed phase diagram, confirming 
that the cubic phase must be a high-temperature ground state of BMT that is metastable 
at room temperature. 
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II. EXPERIMENTAL 



The single crystals of BMT were synthesized according to the procedure described by 
Galasso and Pinto [TTJ . The powder of BMT was prepared according to two-stage synthesis 
described in Ref . [18] and is essentially the same as used in the neutron diffraction study [19] . 

The experiments were performed on X-ray scattering instruments at the ESRF and SLS. 
At the ESRF the diffuse scattering from BMT was measured at the Swiss- Norwegian Beam 
Line (SNBL) and the XMAS beamline, and at the SLS the measurements were made at the 
beamline X06SA. More details of the measurements and the experimental results for BMT 
are given in Ref. [8] . In this paper we concentrate on the unusual diffuse scattering observed 
from BMT and typical results are shown in Fig. [TJ Similar results were obtained from 
the SLS diffractometer X06SA, but from the XMAS data the scattering contained sharp 
peaks at the (H+l/3, K+l/3,L+l/3) positions showing the structure was trigonal. A BMT 
powder sample was also studied at the SLS on the beamline X04SA. It was revealed that 
the trigonal satellites were significantly broader than the intense peaks of the cubic phase 
showing that the domain size was about 20 nm for the trigonal structure. These results 
support the suggestion that our powder and the surfaces of relatively large BMT crystals 
consist to some extent of trigonal-structured domains while the bulk of crystals is cubic and 
shows the diffuse scattering of Fig. [1} 

This paper concentrates on the discussion and modeling of the unusual cube-shaped dis- 
tribution of scattering and on how the different structures of BMT can exist. It is suggested 
that the structures are unable to relax when rapidly cooled from high temperature. 

BaMgi/3Ta2/303 is cubic with lattice parameter a = 0.404 nm. The A cation (Ba) is 
located within the unit cell in (1/2,1/2,1/2), the mixed B-site is in (0,0,0) and the three 
oxygens occupy the (0,1/2,1/2), (1/2,0,1/2), (1/2,1/2,0) positions. 

III. IDEALIZED MODEL OF DIFFUSE INTENSITY 

In this section a simple mathematical model of the diffuse intensity is described. The 
model is simple because each atom is assumed to be a point scatterer , with element Xx 
having a constant scattering length bx x - This means neglecting the atomic form factors and 
Debye- Waller factors, whose effect is just a smooth decay of the intensity at large transferred 
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FIG. 1: (Color online) The diffuse scattering observed for sections of reciprocal space for BMT. 
(Top left) the (100)-(001) section through the origin of the reciprocal space, showing weak diffuse 
streaks (here magnified) of displacive origin; (Top right) the (110)-(001) section through origin 
of the reciprocal space; (Bottom) the (100)-(010) section through (0, 0, 1/2); The square diffuse 
feature, that turn out to be cubes of diffuse scattering, has the intensity largely concentrated in 
the range from (h + 1/3, k + 1/3, 1 + 1/3) to (h + 2/3, k + 2/3, 1 + 2/3). 

momentum values. 
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A. Idealized sharp cubes 



The scattering from BMT has diffuse scattering cubes centered at q = (h, k, I) + (§, §, |) 
with half-side 1/3 (all in rlu). A simple mathematical description of their shape can be given 
as follows. In one dimension, define a box function 

B w {x) = 1 if \x\ < w, 

= otherwise 

(1) 

Then a periodic sequence of such boxes centered on the half-integers is 

®w(x) = B w{x-n-l/2) 

n=— oo 

Multiplying three such terms (one for each coordinate) and adding a constant background 
(that we will show to be necessary) we obtain a model of the X-ray diffuse intensity: 

I dif (q) = I B + I <S> 1/6 {q), (2) 
where 

{1 within the cubes JL- 
= fl E B 1/6 (q a -h a - 1/2) (3) 
OUtSlde a=lh a =-oc 

Here w = 1/6 is the halfwidth of the observed cubes in reciprocal lattice units. Taking the 
Fourier transform of the periodic array of cubes of X-ray diffuse intensity I d if(q) (Eqs. ( 2p )), 
we evaluate the pair correlation function as 



G d if{r) 



f d 3 ge -a,rir., Jd . /(g) 
JR3 



Note that, as Idif is periodic and its Fourier transform is necessarily a weighted set of Dirac 
deltas on the direct lattice [33] 



G dtf (r)= v 2 J2s M S(r-M) 



M 



where if is a positive constant related to the actual scattering lengths (see Sec. Ill C 



Eq. (10)) and Sm are the correlation coefficients. In fact, substituting I^f from Eqs. (2|3) 



we obtain 

sin (7rm a /3) 



G«,{T) = I B 8{r) + £ y E6(r-M)f[ 



27 

M a=l 



-1 



7im a /3 



(4) 



whereas the correlation coefficients are given by 
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^ sin (7rm Q /3) 
nm a /3 



(0,0,0); (5) 



27I B + h , „s 

^(o,o,o) = 27 ^ 2 l b ) 

Here Ai" = (mi,m2,m3) are direct lattice nodes, where the density (and the pair corre- 
lation function) are different from zero. 

The distance r = [m\ + m\ + m|] , so that Gdif(r) ~ r _Qf , with 1 ^ a ^ 3, depending 
on the direction. This is a very rare case of power-law correlation decay of structural 
origin. In quite exotic systems, including spin glasses and nematic crystals, such correlations 
have been described J2U1 EI] and named as quasi-long range order (QLRO). We shall name 
QLRCO (as quasi-long-range cation order) the kind of order exhibited by the BMT sample. 
Two parameters are often used to characterize the degree and extension of the structural 
chemical order. The first is the long-range order parameter of Refs. [22, 23J, see Appendix, 



Sec. VIII C The value of the long-range order parameter X = 1 characterizes periodically 
ordered structures, as the type iii) above, while zero is the value for randomly ordered 
structures, as type i). Conversely, the short-range order parameters of Ref. [21] characterize 
the local degrees of order - supposedly when long-range order is absent. We obtain (see 



Sec. Ill G 1 Sec. IIIG2) that both long- and short-range parameters indicate a substantial 



degree of order, which is an unusual occurrence. 



B. Smoothed cubes 



The model of the diffuse intensity described above can be improved in one significant 
aspect. The cubes as observed are not perfect, as there are intensity fluctuations within 
them. Their edges, especially, are not very sharply defined. This effect can be included by 
convoluting $(g) with a point-spread function (PSF) D(q). Accordingly, the pair correlation 
function Gdif{r) is multiplied by the Fourier transform of the PSF, D(r) as shown in fig. 2 
The latter will be named here the correlation decay function (CDF). This will give a stronger 
decay law for the pair correlation function. This is just what happens with every known 
diffraction phenomenon - as far as infinite crystals do not exist, even the most perfect 
standard crystal will show Bragg peaks convoluted by a suitable PSF and pair correlations 
damped by the corresponding CDF. 
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FIG. 2: (Color online) Effect of a Lorentzian PSF on the cube shape in one dimension. Black - the 
sharp cube profile, blue - the PSF, red - the smoothed cube profile. 

We will not, in the following, burden the notation by including a PSF, but we will 



calculate some of its effects in Sec. |IV| Different possible PSF are discussed in Sec. |VIII A 
Here, as an example, we consider a Lorentzian PSF along the cubic axes (Fig. [2]) 
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D(q) = (nk o y* J] [1 + (qjko) 2 ]' 1 



Its Fourier transform is 



D(r) 



a=l 



n 

Q = l 



-27r|r a |fco 



(7) 



(8) 



The correlation decay length is L = l/(27r/c )- From our data, taking k as the rise of the 
half-width of the cube edges, we can coarsely estimate that L is about 10 unit cells. This 
is well below the coherent domain size, however it remains a considerable distance when 
dealing with diffuse scattering. 

As the PSF D(q) is convoluted with the scattered intensity, the CDF D(r) is multiplied 
by the pair correlation function, that is the Fourier transform of the intensity. Its effect is 
then to limit the spatial extension of the pair correlation function. This is a familiar issue 
in crystallography. Crystal size can be estimated from the broadening of the Bragg peaks, 
due to the limited extension of the crystal - and, in turn, of the pair correlation function. 

We wish to remark here that the 'cube' scattering is quite different from Bragg scattering 
with respect to its response to a limitation of correlation length. In fact, convo luting the 
PSF with a Bragg peak gives the PSF itself. This is the reason for the familiar behaviour of 
Bragg peaks - namely, that their width is inversely proportional to the size of the correlated 
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domain (a good estimate of this is L). The diffuse cubes, however, show a negligible change 
in width, second-order or more in 1/L. This scattering form is therefore fundamentally 
different from Bragg scattering. 

C. B-site order and diffuse cubes 

To explain such diffuse scattering in terms of B-site cation order we split the scattering 
density p{r) as 

p{r) = p A (r) + Ap(r), (9) 

with pa the periodic average and Ap a non-periodic zero-average fluctuation. If the B site 
is the only one that varies, all of Ba and O partial densities are included in p^. As already 
mentioned, we simplify the picture by considering point-like atoms. Then Ap(r) is zero 
everywhere except on the B sites, and for B=Ta Ap(r) = 77; for B=Mg Ap(r) = —27/. If 
&M g and &Ta are the scattering lengths of Mg and Ta, respectively, then 

V = (&Ta - &Mg)/3. (10) 

As Ap assumes only two values (r/ on Ta sites and — 2rj on Mg sites) we can consider 
(see PQ) the site concentration of one of the two species, say Mg, in each lattice node. Clearly 
this can be only or 1. We indicate it with cm, that is 1 if site M is occupied by Mg and 
otherwise. Similarly we take c' M as the Ta local concentration, that is 1 if site M occupied 
by Ta and otherwise. As we have no vacancies, c' M = 1 — cm- Averaging over all sites, we 
have 

(c) = 1/3 and (c') = 2/3. 
Furthermore, as c and d can be only or 1, 

(c 2 ) = ( c ) = 1/3 and (c ,2 > = (c'> = 2/3, 
so the r.m.s. fluctuation of the concentration is given by 

5c = vV) - (c) 2 = V2/3. 
The fluctuation of the density is given by 

Ap(r) = V J2(~ 2c m + c' M )5(r -M)= V ^(1 - 3c M )5(r - M) (11) 

M M 
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D. Correlation function 



In order to study the diffuse scattering from BMT, we need to evaluate the pair correlation 
function of Ap: 



G dif (r) = Q 



\f d 3 r'Ap(r')Ap(r' + r) = ri 2 \2s M S(r-M) (12) 



Here 

n = Nv c = N 

is the crystal volume, N the number of lattice nodes within it, v c the unit cell volume 
(taken as 1 hereafter, hence the last equality). We shall equate this to the pair correlation 
function derived from the data, see Eq. Q. The rightmost form is because Ap is a weighted 
Dirac lattice, and this property is invariant upon self-convolution. This is the same form as 
Eq. Q, so clearly now we have to calculate the Sm pair correlation coefficients for a given 



fluctuation density (Eq. (11)) and to compare them with those obtained from the Fourier 



transform of the data (Eqs. ( 5p )). 



Substituting Eq. (11) into Eq. (12) gives 

,2 



G*f(r) = ^^^(l-3c i )(l-3c jL+M )<5(r-M) (13) 



M L 



where, comparing with Eq. (12), we obtain 



9 \ -< 

Sm = -1 + -jy y J c l c l +m = -1 + 9 (c L c L+M ) L ; ^(ooo) = -1 + 9(c 2 ) = 2. (14) 

L 

Now, the probability of finding a pair of Mg atoms separated by a lattice vector M is 
given by 

P Mg -Mg(M") = (c L c L+M ) L = + 9 M (15) 

Similarly we can calculate the probabilities of all the other pairs (Mg-Ta/Ta-Mg and Ta-Ta) 
for any spacing vector M, and the result depends only on the Sm correlation coefficients: 

PTa-Ta(M) = <(1 - C L ){1 - C L+M )) L = ^ + (c L C L+M ) L = ^ + ^ M 

P mix (M) = (c L (l - c L+M )) L + ((1 - c L )c L+M ) L = g - 2 (c L c L+M ) L = ^— ^16) 

where the subscript 'mix' refers to Mg-Ta pairs, in any order. Therefore, if we know the cor- 
relation coefficients we can calculate all the properties that depend on the pair interactions. 

13 



E. Scale and background 



Comparing Eq. Q with Eqs. ( 13|14 ), considering r = 0, we obtain 



S, 



(000) 



Ib 



27 



54r] 2 - nil 



We define a parameter t, so that if ^ t ^ 1, then 

I B = 2rf{l-t)- h 



hArft 



(17) 



(18) 



This result describes how much of the total diffuse intensity in one Brillouin zone is a 
constant background and how much comes to form the cubes of diffuse scattering. 



As described by Eq. (15), (1 + 5m)/9 is a probability, hence it must range between and 



1. In order for it to be positive for all M, 



Sm iS — 1. 



For M = (1, 0, 0), this constraint - using Eq. ^ - implies 



t < 



7T 



3^ 



0.6046, 



(19) 



(20) 



which tells us that a consistent fraction - at least about 40% - of the diffuse intensity must 
be spread as a constant background while the rest forms the diffuse cubes. This is also a 
lower bound only, so possibly the background fraction is even larger. 



F. Electrostatic energy as a function of the pair correlation function 

we have calculated the probability of finding pairs of atoms of equal or 
different species (limited to the B sublattice) for each separation vector M. That is all we 
need to calculate the electrostatic energy. In fact, the structure consists of the fixed part 
(the A sublattice and the oxygen sublattices) and the variable part (the B sublattice). A 
result of Ref. [25] is that the electrostatic energy does not depend on the interaction terms 
between the different sublattices; the fixed part then has no influence except for a constant 
term and one needs to calculate only the interactions between atoms of the B sublattice. 
Similarly to the scattering density (Eq. (|9|), we can separate the charge density p e (r) into 
the the average part and fluctuating part 

Pe(r) =Pe,A( r ) + A Pe(r), (21) 
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In Eqs. (15|16 



with p e A the average component and Ap e a non-periodic zero-average fluctuation. Again, 
on Ta-occupied sites we will have Ap e = r\z and on Mg sites, Ap e = —2r] Z with 



Vz = (Z Ta - Z Mg )/3. (22) 

If we neglect covalent bonding and assume the formal charges = +5, Zyi g = +2, then 
we have rjz = 1. It can be shown [23] that only the fluctuation density Ap e - that is the 
only term that varies between the different possible phases of BMT - needs be included in 
the calculations. Therefore, the electrostatic energy is simply expressed - using Eqs. ( 15|16 ) 
- as 

E = ^Vz_ 4P Mg -M g (M) + P Ta -Ta(M) - 2P mix (M) = A Sm_ 

8iree ^ \M\ 87ree ^ \M\ [ ' 

the prime meaning that M = (0, 0, 0) is excluded from the summation. An additional factor 
of 1/2 was included in order not to count twice each bond. 



G. Monte-Carlo solution 

At this point we know that a distribution of B-site cations whose correlation coefficients 
are given by Eqs. (5|6) would produce a diffuse intensity as described by Eqs. (2|2), that well 



matches the observed diffuse scattering cubes. We cannot, however, prove analytically that 
a distribution of two different cations on a cubic lattice with respective abundances of 1/3 
and 2/3 exists that can satisfy Eqs. (5|6) (or alternatively, produce a diffuse intensity as in 



Eqs. (2p)). In order to prove that such a distribution exists, we decided to try to fabricate 
one or more by a Monte-Carlo simulation. In this way we could also get some structural 
insight by analyzing the configurations so obtained. 

A large cube (27 x 27 x 27 unit cells with the total number of atomic sites iV = 19683) 
was placed on the lattice nodes with a 'Mg' atom - represented by a scattering length —2 - 
or by a 'Ta' atom - represented by a scattering length 1. These values are chosen so that 



the average density (cf. Sec. Ill C Eq. ([9])) is zero and only the fluctuation density Ap 
is not zero; moreover, we set r] — 1 for simplicity. The FFT power spectrum was evaluated 
from the intensity Idif(q) in one Brillouin zone. The atoms were moved around until a good 



agreement with Eqs. (2j3) was achieved. The simulation was very successful. In Fig. [3] a) 



a section of the calculated diffuse intensity is shown for one of the optimal configurations 
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found. The result of the calculations compares well with the experimental data. (Eqs. 
Other sections (not plotted) are consistent with the one shown. 




FIG. 3: (Color online) a): The diffuse intensity calculated in one Brillouin zone from the Monte 
Carlo solution. The result shown in the (h,k, 14/27) plane. 

b): Red bars: correlations Sm - f° r M = (A;, 0,0) - from the Monte Carlo calculation. As allowed 
by the cubic symmetry, values have been averaged over M = (±fc,0,0), (0, ±A;, 0), (0,0, ±fc), and 
plotted vs. k. Blue line: correlations Sm calculated (on the same M points) using Eqs. ( pp ). 
This is plotted as a double envelope curve (because of the factor — l fc ) to guide the eye. Green line: 
same as the former, multiplied by a CDF exp(— k/L) as from Eq. ([8]) with L = l/(27rfco) = 10. 



1. Pair correlations and long-range order parameter 

Several equivalent configurations were constructed as above and used to evaluate the 
pair correlations. Firstly, the ratio Iq/Ib for the solutions was found, and it resulted that 
Iq/Ib ~ 3.0. This gives a value of the intensity partition parameter t ~ 0.1 - actually much 
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lower than its upper bound of 0.6. From that, exploiting Eq. (29) in the Appendix, the 
value of the long-range order parameter was extracted asX~ 0.407. We can then evaluate 
I B = 1.8r/ 2 , Jo/27 = 0.2 V 2 . 

Fig. [3|d) shows the correlations Sm from the Monte Carlo compared with those evaluated 
above, for M along the family of (1,0,0) directions). The good agreement with values 
calculated from Eqs. (5|6), with and/or without an additional PSF with L = 10, shows more 
quantitatively that the Monte Carlo converged simulation converged towards the observed 
pattern. 



2. Crystal chemical information and short-range order 

In order to display some crystal-chemical information in the MonteCarlo-obtained QL- 
RCO configurations, we directly evaluated the probabilities that a Mg (respectively, Ta) 
atom has a Mg (respectively, Ta) neighbor at a distance k along any of the equivalent 
(1, 0, 0) directions. The results are shown in Fig. [4] and have been calculated as follows. 
Each site in a cluster is treated as the central atom. The fraction of atoms of a given species 
over the neighbors having distance k along the family of (1, 0, 0) directions was calculated. 
Periodic boundary conditions were used. The fractions for all four chemical pairs where the 
atom in the origin is Mg or Ta and the neighbor is Mg or Ta, respectively. In Fig. |4| we plot 
the Ta-Ta, Ta-Mg, Mg-Ta, Mg-Mg (kOO) chemical pair fractions for the QLRCO cluster. 
For comparison we add the values for a cluster with randomly assigned Mg or Ta at each 
site, fulfilling a composition Mg 1 / 3 Ta 2 /3. These of course do not depend on the choice of 
the atom at the origin, and amount to 1/3 when the neighbor is Mg and 2/3 when it is Ta. 
From the values of the Ta-Mg (or Mg-Ta) chemical pair fractions we immediately obtain 
the value of the short-range order parameter for the first coordination shell [24, 26jfor the 
QLRCO phase. The parameter a(l) is found to be ~ —0.21. For this composition it has 
a maximum (in absolute value) at —0.5, for the fully ordered state, while it is zero for the 
random state. Combined with similar results obtained for the long-range order parameter 



in Sec. HI G 2 this gives a further justification for the name 'quasi-long range order' 
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(Ta,Mg)-Mg random 
(Ta,Mg)-Ta random 
Mg-Mg QLRCO 
Ta-Mg QLRCO 
Mg-Ta QLRCO 
Ta-Ta QLRCO 



FIG. 4: (Color online) Blue: Mg neighbor fractions for a central Mg atom in the QLRCO 
Monte-Carlo-generated cluster, averaged along all (kOO) separation vectors. These correspond 
to 3P Mg -Mg((fc00)) (c/. Eqs. ( |15|16[ ), the factor 3 weighs the probability of the central atom 



being Mg). For the same cluster: red: Mg neighbor fractions for a central Ta atom, equal to 
(3/4)P m i x ((A;00)) (here factor 3/4 is obtained as product of 1/2 - because P m ix includes both bond 
directions, and 3/2 - inverse of the probability of the central atom being Ta). Green: Ta neighbor 
fractions for a central Mg atom, (3/2)P m ; x ((A;00)). Brown: Ta neighbor fractions for a central Ta 
atom, (3/2)Pxa-Ta((^00)). For a cluster with randomly assigned Mg or Ta at each site, but with 
the exact composition Mg 1 / 3 Ta2/3, we pl°t i n black the Mg neighbor concentrations of any central 
atom; in violet, the Ta neighbor concentrations of any central atom. These values can be calculated 
as above, simply substituting S^oo) = 0- 

IV. ELECTROSTATIC ENERGY CALCULATION RESULTS 

In this section and the next we shall demonstrate how the electrostatic energy and the 
configuration entropy, as calculated from the structure and pair correlations for the different 
possible BMT phases is sufficient to explain the phase diagram of BMT. In particular, the 
fact that the QLRCO phase is metastable at room temperature, and the fact that the NaCl 
phase, although common in similar perovskites, has never been found for BMT. 

The electrostatic energy, is not the only component of the crystal energy, but has been 
previously identified as the driving force of the cation ordering in complex perovskites. 
At least it has been shown [25] that electrostatic energy-based predictions describe quite 
effectively the ordering behaviour of a large class of perovskites. The electrostatic energy 



was evaluated in a similar fashion as in Ref. [25J. As already mentioned in Sec. Ill F only 
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the fluctuation term relevant to the B-cation sublattice need to be calculated, because the 
other contributions can be considered either as fixed or as zero. 

We have used a relative permittivity e = 1 and a cubic cell side a = 4 A and the formal 
charges Ta +5 , Mg +2 , Ba +2 , 0~ 2 which gives r/z = 1. It is easy to scale the results to more 
realistic values of the parameters later on. 

The electrostatic energy was calculated using different ways. In all cases, the zero of 
the energy was fixed to the energy of a reference perovskite - BaTiOs - that has a +4- 
charged cation (equal to the average charge of Mg 1 / 3 Ta 2 /3) on the B site. Electrostatic 
energies were calculated - apart from the reference BaTiOa perovskite, labeled uniform, 
for all the different B-site ordering discussed in this paper. The BMT phase (we label it 
'QLRCO' for conciseness), the structure with NaCl-type partial B-site order (label 'NaCl'), 
the trigonal phase with full B-site order (label 'P3ml'), and the phase having fully random 
B-site occupancies (label 'Random') have been calculated. 

As a first method of calculation, we created large clusters of atoms and computed the 
interactions directly, summing up all atom pair contributions. Calculations of this kind 
were performed for clusters of increasing size up to 27 x 27 x 27 cubic cells, at which 
limit the results were well stabilized and appeared to have reached their asymptotic value 
(with residual oscillations smaller than 0.1 eV/formula unit). Fractional O- and Ba- site 
occupancies at the boundary have been used to keep the neutrality of the clusters. This 
fact |27J, together with the cubic cluster symmetry, explains the very rapid convergence of 
the electrostatic energy sum. 

As a second approach, the Madelung method [27] was used. This gave consistent results 
within < 0.1 eV/formula unit. 



Thirdly, we use the knowledge of the correlation coefficients Sm (cf. Sec. HID). In this 



case, the electrostatic energy is simply expressed as in Eq. (23). Such sums were evaluated 
in all cases from the analytical Sm- Again, we evaluated them on cubes of increasing size 
(up to 500 x 500 x 500 unit cells) in order to verify the convergence - which was always 
excellent. Calculation times in this case were also extremely fast - less than one minute 
on a single PC for the largest cubes with « 10 8 nodes. With this latter method for the 
QLRCO case was also evaluated - the effect of a PSF spreading the diffuse scattering cubes 
(see Fig. [5]). It turns out that even with small correlation lengths (L ~ 10 unit cells) the 
energy was within few % from the value as L — > oo. 
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FIG. 5: (Color online) Electrostatic energy - Eq. (23) - for the QLRCO structure, when correlations 
are reduced by various decay functions (cf. Sec. Ill, Sec. VIII A). 



For QLRCO, NaCl and Random cases (which are not deterministic) we have repeated 
the calculation for many configurations, evaluating the average values and their standard 
deviations. In all cases the number of trials was chosen high enough for the standard 
deviations to be below 0.1 eV/formula unit. Permittivity was taken as e — 1, and the 'true' 
energy levels can be obtained by dividing the results by the permittivity of the material. 
The results for the uniform structure could be compared with those evaluated using the 
Madelung constants (see (27]) and turn out to be quite accurate. Also the electrostatic 
energy of the 'P3ml' structure - calculated by the Madelung method by assuming that 
the site geometry remains cubic is found in Ref. [2H] and compares well with our results. 
Furthermore we are not interested in absolute values but in comparative ones, and these are 
numerically much easier to compute than the absolute ones. For each structure we evaluate 
AE s t rU cture = E 'structure — -^BaTicv The results show that, with respect to the Uniform 
structure, in order of increasing stability: 

• the Random structure has the same energy within 0.1%, lAi^Randoml ^ 0.1 eV/formula 
unit; 

• the QLRCO structure has 1.1% lower energy or about APqlrco = —1-902 eV/formula 
unit; 

• the NaCl structure has 1.8% lower energy or about AE Na _ci — —3.229 eV/formula unit; 

• the P3ml structure has 3.0% lower energy or about AE P ^ ml = —5.295 eV/formula 
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unit. 



Noteworthy is the level spacing ratio 



c 



A-EqlRCO — Ai^NaCl 

A£ NaC i - AEpz ml 



0.642 < 1 



which is smaller than 




P3ml 



'NaCl 



1.56 > 1 



There are several factors which will influence the absolute values of these energy levels, 
and the relative values of ( and £' as well. These are: 

• The energies have to be divided by the material's relative (static) permittivity e, that 
can be large in perovskites. For BaMg 1 / 3 Ta 2 / 3 3 it is found that e ~ 24 |2j [29]. This 
factor would not change the ratio. 

• The electrostatic energy is evaluated using the full formal charges of the ions. Cova- 
lency corrections may play a role. A value of r/z = 0.913 can be extrapolated from 
the valence sum corrections evaluated for BMT in Ref. [2J. This entails an absolute 
change of the energy levels, not affecting the ratio (. 

• The electrostatic energy is inversely proportional to the cubic lattice parameter a. 
Appropriate scaling is needed. However, for the known BaMg 1 / 3 Ta 2 / 3 3 phases, P3ml 
and QLRCO, a (evaluated from the cubic main reflections) is the same within one part 
in 10 3 [2]. This effect cannot change the ratios. 

• Geometric distortions are present in the ordered NaCl and P3ml structures, but they 
affect the distances on few parts in 10 4 [2J, and as their effect on the energies is of 
comparable order, we can neglect them altogether. 

• A PSF is observed, leading to smoothed cube edges. This term is possibly 'accidental', 
depending on the detailed crystallization kinetics, so it will not be considered quan- 
titatively in the following. Its effect is also small enough (see Fig. [5]) that it cannot 
qualitatively affect our results. 
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V. ENTROPY 

The crystal structure can be stabilized by the entropy For disordered structures config- 
urational entropy plays an important role. A mixed-occupancy atomic site with occupancy 
x for one species and 1 — x for another gives an entropy contribution of 



In the P3ml structure no mixed sites exist, therefore they have no entropic contribution to 
the Gibbs free energy. There will be an entropy term related to the domain distribution, 
but that is smaller by many orders of magnitude. 

The maximal configurational entropy per formula unit is expected in the random structure 
and it amounts to 



For the QLRCO structure, as discussed, the entropy value is slightly smaller, but the 
difference can be evaluated numerically (in a simple but tedious way) from the correlation 
coefficients and the result is that the correction is negligible, so the entropy is ~ So- 

For the NaCl structure, only 50% of the sites are mixed so the entropy is So/2. 

One more important contribution is the phonon entropy. A different B-site degree of order 
can in principle affect the phonon entropy. However, this effect is much smaller than the 
configurational term. This is supported by the scarce literature, see e.g. [30] • The specific 
heat C v for BMT has been measured [31], however the phase was not fully characterized 
at the time (being, however, either QLRCO or P3ml). We can only compare the entropy 
evaluated from the experimental specific heat with data for Pb containing relaxor perovskite 
(PbMg 1/3 Ta 2 /303-PMT) [3H [32]. The phonon entropy of PMT extrapolated to T> 6^ 
exceeds that of BMT by 59 /zeV/K. This large difference can be explained either by the 
polar nano-regions which are only present in PMT [3T] or by the difference in the masses of 
Pb and Ba. 

VI. DISCUSSION AND CONCLUSIONS 

Entropy is important at high temperature for stabilizing different structures. For instance, 
considering the various free energies as calculated in the previous section, we can construct 



S 



k B [x log (x) + (1 — x) log (1 — x)] 



(24) 




(25) 
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FIG. 6: Schematic phase diagram of BMT. 

a phase diagram of BMT (Fig. [6]). We have taken, as previously stated, realistic values of 
e = 24 (assumed constant with T), r]z = 0.913, and a relative linear expansion coefficient 
of 9xl0 -6 . We have not included any contribution from imperfect ordering (namely, the 
effect of a CDF limiting the spatial extent of correlations) because it is not easy to evaluate 
quantitatively and it is not an intrinsic effect but rather a kinetic effect while its contribution 
is small (Fig. [5]) and the leading approximation would only compress the energy range. 

The entropy term is canceled by the energy difference between QLRCO and P3ml phases 
at about 1600 °C, that is about the melting point. Above this temperature the QLRCO state 
would be stable while the P3ml state would be stable below. We can also see from Fig. [6] 
that there is no region that is stable for the NaCl-ordered phase, and in fact, this phase has 
never been observed in BMT, at variance with PMT [21 H]. 

Due to the slow kinetics associated with cation ordering, chemical ordering is usually 
achieved only after annealing for a long time at very high T, where the equilibrium state 
may be dictated by the entropy. Cooling - due to the slow kinetics - cannot change the high-T 
equilibrium state so that it would remain preserved (metastable) at room temperature. 

The random structure is always less stable than the QLRCO structure, because of the 
tiny entropy difference and the large energy difference that favours the latter structure. 
Therefore the random structure is a disordered metastable state. The ground state - or the 
thermodynamic equilibrium state at K - is the P3ml structure. The true high-temperature 
thermodynamic equilibrium state is QLRCO, because the electrostatic energy difference with 
the P3ml structure is overcome by the entropy term at high T, even if the dielectric constant 
is of moderate size. The NaCl structure is intermediate between these two. Given that its 
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Gibbs free energy decreases with increasing temperature half as fast as the free energy of 
the QLRCO state, the ratio £ becomes crucial to determine if the NaCl order can be the 
thermodynamic equilibrium state in some temperature range. In fact for ( < 1 the Gibbs 
free energy of QLRCO is always lower than that of NaCl structure, while for £ > 1 there 
exists a temperature interval where NaCl structure is the equilibrium state. The QLRCO 
structure is favoured above that temperature and the P3ml structure below. 
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VIII. APPENDIX 



A. More choices of point-spread functions 



A particular type of PSF/CDF pair has been shown in Eqs. (7p). In this section some 



different PSF-CDF pairs will be discussed. PSF's are normalized so that their sum is 1. 
Firstly, as in Eqs. (7l8), we consider the symmetry of the PSF to be cubic. A Lorentzian 



PSF and the corresponding exponential CDF are presented in Eqs. (7l8). A simpler PSF 



would be a sine function, corresponding to a box function CDF. The negative parts of the 
sine tails are, however, a problem. Therefore we consider a sine squared function: 

PSF: D(q) = (nk y 3 l[ S1 f ^ ^ ; CDF: D(r) = - n\r a \k )hox(n\r a \k ) 

a=l lW«oj a=1 

(26) 

Correspondingly, if the symmetry is spherical, the PSF/CDF are given by the radial 
Lorentzian or exponential: 

PSF: D(q) = 7i- 2 k 3 [l + (q/k ) 2 Y 2 ; CDF: D(r) = e~^ rk ° , (27) 
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and the analog of the squared sine/box pair: 



PSF: D(q) = 3(2n 2 k% 



2 7„3\-l 



sin (q/kp) - (q/kp) cos (q/kp) 
(g/fco) 3 



CDF: D(r) = (l+nrko/2)(l-nrkp) 2 box(nrk ) 

(28) 



B. Correlation coefficients for ordered and semiordered BMT phases 

In the case that the B-site is fully disordered, i.e. when all c& are randomly assigned 
5(000) = 2, but 

S M = -1 + 9(c) 2 = whenever M ^ (0, 0, 0) 

so the Fourier transform I(q) is a constant, I(q) = 2rf. 

The trigonal superstructure has as the 3-fold axis one of the equivalent (1, 1, 1) cubic 
axes. This means that it exists in four different orientations, where the 3-fold axis T is 
T = (1, 1, 1), T = (1, -1, -1), T = (-1, 1, -1), T = (-1, -1, 1), respectively. The B-site 
superstructure is obtained by positioning in a lattice node M on a Mg atom if M • T = 
modulo 3. Therefore, we will have pure Mg planes alternating with pairs of Ta planes, 
stacked orthogonally to T. In any crystallite, domains of the four different orientations will 
be present in equal volumes, therefore conserving the global cubic symmetry, in agreement 
with Landau's theory. The correlations can be calculated for each orientation and then 
averaged. However, one can also calculate them from the scattering pattern, as an average 
over domains. The superstructure satellites - considering equipopulated domains - are eight 
equal Bragg peaks (each with intensity ip3 ml ) at (±1/3, ±1/3, ±1/3). Then 

S M = (_i)*+i2*^pi 
V 

where k is the number of components of M that are equal to zero, or k = #{m a = 
modulo 3, a = 1,2,3}. From the condition S'(ooo) — 2, we obtain ipg ml = — , so that 

S M = (_i)fc+i 2 fe - 2 

The hypothetical NaCl-type superstructure is simpler. There is only one Bragg peak at 
(1/2,1/2,1/2) with intensity iNaCi phis a small constant background I b due to residual ran- 
domness on half of the sites. So one obtains 



-, ,„ 7 ,„ . ,„ :! I\nC . _ I\m + h 
>m — y— -U — 2~, *J(ooo) — 2 
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As half the sites contribute to the background and the occupancy of the random sites is 
always divided as 1/3 - 2/3, lb = i] 2 and again using the condition S(ooo) = 2, we have 

^Naci = V 2 - 

C. Long-range order parameter 

In this section we show how the former results can be related to the formalism of the long- 
range order parameter. A simple way of describing systems where there is a fully ordered 
state is by the long-range order parameter [22l 123] . Here we will not go to all the detail 
as it is a well-known theory. For BMT the ordered state is the trigonal modification. Here 
the cubic B-site lattice is divided into 3 sublattices, denoted (a),(b),(c) for convenience of 
which one is occupied by Mg and two by Ta. There are four ways of choosing the three-fold 
axis of the trigonal structure and for each three permutations of the Mg sublattice between 
(a),(b),(c). Suppose that we chose one axis and that the Mg sublattice is (a). In the ordered 
state all Mg atoms are actually on the (a) sublattice; if we have N total atoms, we will 
have A/3 Mg atoms in the (a) sublattice. In the random state, the Mg concentration on 
any site will be the same, so we will find only iV/9 Mg atoms in the (a) sublattice. We 
can parametrize this situation if we say that the number of Mg atoms on the (a) sublattice 
is given by A(l + 2A)/9 where A is a parameter that is 1 in the ordered state, in the 
random one and assumes intermediate values for all intermediate states. Using the same 
A, we can parametrize the number of Ta atoms on the (b) + (c) double sublattice, and all 
cross-frequencies. More involved calculations (see [23] for details of the method) lead to the 
expressions for the nearest-neighbor pairs frequencies as a function of A. If the interaction 
is fundamentally short-ranged - or if the interactions beyond the first neighbors are not 
very important - one can then use different interaction hypotheses and calculate the internal 
energy. Finally, from the stochastic occupancy of each site, it is possible to evaluate the 
configurational entropy. The nearest-neighbor bond fractions turn out to be 

ui = -(1 — A 2 ) Mg-Mg pairs 
9 

v 2 = -(4 — A 2 ) Ta-Ta pairs 
2 

zy 3 = _(2 + A 2 ) Ta-Mg pairs 
9 
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As nearest-neighbor bond vectors are just the (1, 0, 0) family, it is easy to see that 

U l = TP 2^ C L C L+(l,0,0) ~ 



Similar relations for the other z/s can be given using the c' L = 1 — cl's. Since X 2 = — 5(1,0,0) > 
and using Eqs. ( 5|18 ), a relationship is found between X and the background parameter t 



(see Sec. HIE) 



X 



71 



1/2 



(29) 



The correlation coefficients are evaluated as in Sec. |VIIIB[ while X = 1 for the trigonal 
phase but also for the partially ordered hypothetical NaCl-type phase. Obviously X = 
for the random B-site phase, while for the QLRCO phase, using the t ~ 0.1 value estimated 
from the MonteCarlo runs, we have X = 0.407, which is well above the value of zero for 
the completely random structure. This - combined with similar results obtained for the 



short-range order parameter of [21] in Sec. IIIG2 - is a further justification for the name 
'quasi-long range order'. 
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